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ABSTRACT 

The spatial morphology, spectral characteristics, and time variability of ultracompact H II regions 
provide strong constraints on the process of massive star formation. We have performed simulations 
of the gravitational collapse of rotating molecular cloud cores, including treatments of the propagation 
of ionizing and non-ionizing radiation. We here present synthetic radio continuum observations of H II 
regions from our collapse simulations, to investigate how well they agree with observation, and what 
we can learn about how massive star formation proceeds. We find that intermittent shielding by dense 
filaments in the gravitationally unstable accretion flow around the massive star leads to highly variable 
H II regions that do not grow monotonically, but rather flicker, growing and shrinking repeatedly. 
This behavior appears able to resolve the well-known lifetime problem. We find that multiple ionizing 
sources generally form, resulting in groups of ultracompact H II regions, consistent with observations. 
We confirm that our model reproduces the qualitative H II region morphologies found in surveys, with 
generally consistent relative frequencies. We also find that simulated spectral energy distributions 
(SEDs) from our model are consistent with the range of observed H II region SEDs, including both 
regions showing a normal transition from optically thick to optically thin emission, and those with 
intermediate spectral slopes. In our models, anomalous slopes are solely produced by inhomogeneities 
in the H II region, with no contribution from dust emission at millimeter or submillimeter wavelengths. 
We conclude that many observed characteristics of ultracompact H II regions appear consistent with 
massive star formation in fast, gravitationally unstable, accretion flows. 



1. INTRODUCTION 



Ultracompact (UC) H II regions ha ve radii R < 0.1 p c 
and high radio surface brightness ( |Churchwell| [2002 ) . 
They have a characteristic distribution of mo rphologies 
( |Wood fc Churchwell|1989| |Kurtz et al."||1994p , and usu- 
ally are found associated with one another and with 
compact H 1 1 regions (e.g. 



Welch et al. 1987; Gaume 



Claussen|1990 Mehringer et al.|1993[|Kim & Koo 2001). 



Lheir spectral energy distributions (SEDs) can reflect a 
transition from optically thin to optically thick emission, 
but often sh ow an anomalous intermediate w avelength 

Recent 



dependence (Franco et al. 2000 



Lizano 



2008) 



radio continuum observations have suggested that ultra- 
compact H II regions can contract, change in shape, or ex- 
and anisotropically over intervals of a s little as ~ 10 yr 
Franco-Hernandez & Rodriguez 2 004| |Rodnguez et al. 
U(J7| |Galvan-Madrid et al.||2008p . 



The observed brightness and size of UC H II regions 
require that they be ionized by massive stars of spectral 
type earlier than B3. If the regions were to expand at 
the sound speed of ionized gas, q ~ 10 km/s, they would 
have lifetimes of roughly 10 4 yr. Less than 1% of an OB 
star's lifetime of a few million years should therefore be 
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spent within such a region, so the same fraction of OB 
stars should now lie within them. However, surveys find 
numbers in our Galaxy consistent with over 10% of OB 



stars being surrounded b y them ( |Wood fc Churchwell 
1989; |De Pree et al.||2005[), or equivalently, lifetimes ol 



~ 10 yr if this model is correct. 

A number of explanations have been proposed for this 
lifetime problem, in cluding confinement in cloud cores 
by thermal pres sure ( De Pree et al.fl9 95; Garcia- S egura | 
Franco|[l996l) or turbulent pressure IXie et al.||r996), 



ram pressure con finement by infall ([Yorke 1986; Hollen- 



bach et al.||1994|) or bow shocks dVan Bur e n. et aL| 1990 ; 

Hoare 2006), champagne 



Mac Low et a ETT991 all Arthur 

flows (Bodenheimer et al. 1979[ |Garcia-Segura fc Franco 
1996 ; Arthur fc Hoare 20(J6|), disk evaporation ( [Hollen- 
bach et~a l. 1994), "and mass-loaded stellar winds (IDys'on 



et al.||1995l IKedman et~aT1[T996l |Williams et al.||1996| 
Lizano et al. ||1996|), but most have been argued to have 



major flaws flVlac Low et al.||2007 ). 



We have modeled accretion on to an ionizing source 



using three-dimensional simulations (Peters et al. 2010 
hereafter Paper I). These calculations suggest that accre- 
tion can indeed explain the l i fetime problem, but in an 
unexpected way. |Keto| ( |2002| |2007| ) has argued that ul- 
tracompact and hyp er compact H II regions are simply the 
ionized portion of an accretion flow. However, massive 
stars require accretion at rates exceeding 10 -4 M yr -1 
flBeuther et al.||2002| jBeltran et aT]|2006| ) to reach their 
final masses be fore exhausting their nuclear fuel ( |Keto| 
|fc Wood||2006[ ). The result is gravitational instability 
during collapse, leading to the formation of dense gas fil- 
aments in the rotating, collapsing flow, along with dozens 
of accompanying stars (Paper I). The strongest sources 
of ionizing radiation orbit through the dense filaments 



2 



repeatedly, accreting mass efficiently when they do. The 
filaments absorb the ionizing radiation locally, though, 
when this happens, shielding the rest of the H II region 
for long enough for it to recombine. As a result, the size 
of the observed H II region remains independent of the 
age of the ionizing star until the surrounding secondary 
star formation cuts off accretion on to the primary and 
a compact H II region begins to grow around it. 

In this paper we consider in more detail than in Pa- 
per I whether our models of ionization interacting with a 
gravitationally unstable accretion flow can reproduce the 
observations of ultracompact H II regions summarized at 
the beginning of this section. We show how H II re- 
gions fluctuate in size as their central stars pass through 
density fluctuations, and demonstrate that our models 
qualitatively reproduce all morphologies observed for ul- 
tracompact H II regions, even giving general quantitative 
agreement with the distribution of differe nt mor pholo 



gies observe d by Wood & Churchwell ( 1989 ) and |Kurtz 
et al. (1994). Our models also otter a natural explanation 
tor the observed clustering of ultracompact H II regions. 
We further demonstrate that they reproduce observed 
SEDs, and provide natural explanations for the anoma- 
l ous SEDs observed for some ultracompact H H re gions 
flLizano||20Q 8| |Beuther et aT|2004| |Keto et al ||2QQ8| ) . 

In Sect. [2] we describe our methods lor modeling ultra- 
compact H II regions and simulating observations, while 
in Sect. [3] we describe the results of our work relevant for 
this paper. Finally, in Sect. [4] we draw conclusions. 

2. NUMERICAL METHOD 
2.1. Simulations 

We present three-dimensional, gas dynamic, simula- 
tions with radiation feedback from ionizing and non- 
ionizi ng radiation. We u se the FLASH adaptive-mesh 
code ( |Fryxell et al. |2QQQ ), modifie d to include a hybrid- 
charactenstics raytracmg method (Rij khorst et al.| 2QQ6) 
to solve the radiative transfer proble m. The proto stars 
are modeled by sink particles ( Federrath et al.|2010 ) that 
are coupled to the radiation module via a protostellar 
model (Paper I). 

We simulate the collapse of a massive core with a mass 
of 1000 M®. The core has constant density within the 
sphere with r < 0.5 pc, while further out the density falls 
off as r -3 / 2 . The gas temperature is initially T = 30 K. 
The core begins in solid body rotation with a ratio of 
rotational to gravitational energy j3 = 0.05. 

We use an adaptive mesh with a cell size at the highest 
refinement level of 98 AU. Sink particles are inserted at 
a cut-off density of p cr i t = 7 x 10 -16 gcm -3 and accrete 
all gas above p cr i t within an accretion radius of r s i n k = 
590 AU if it is gravitationally bound to the particle. 

We analyze two simulations. In the first simulation 
(run A), we only follow the evolution of a single star and 
suppress the formation of any secondary stars. We do 
this by introducing a dynamical temperature floor 



Gfi 
7rk B 



p(nAx) 2 



(i) 



with Newton's constant G, mean molecular weight /1, 
Boltzmann's constant /cb, local gas density p, and cell 
size Ax. This temperature floor guarantees that the 
Jeans length is always resolved with n cells. We must 



set n> 4 to prevent artificial fragmentation ( |Truelove 
et al.|199"7 ). In the second simulation (run B), we do not 
apply the temperature floor, instead allowing the forma- 
tion of secondary sink particles. 

More details on the simulation method as well as a de- 
tailed description of the simulation results can be found 
in Paper I. 

2.2. Generation of Free-Free Emission Maps 

Radio continuum emission from H II regions around 
massive stars at wavelengths A > 0.3 cm (y < 10 11 Hz) is 
predomin antly caused by free-free emis sion from ionized 
hydrogen (Gordon & Sorochenko 2002) . Since scattering 
can be neglected lor this problem (p^raus| 1966} |Gor don fc 



|Sorochenko|2002| ) , the equation 01 radiative transfer can 
be readily integrated. First, we calculate the free- free 
absorption coefficient of atomic hydrogen 



0.212 



( ne V 

VI cm -3 / 



1 K 



-1.35 



(l Hz) 



-2.1 



cm 



( 2 ) 

with number density of free electrons n e , electron tem- 
perature T e an d frequency v. Since the electrons ther- 
malize quickly (Dyson & Williams 1980[) , we can take 
the gas temperature T — l e . Given the absorption coef- 
ficient, the optical depth at distance r from the edge of 
the domain is then 



a v ds. 



(3) 



The radiative transfer equation in the Rayleigh- Jeans 
limit then leads to the brightness temperature 



T{r v ) 



Jo 



(4) 



The resulting map of brightness temperatures can be 
converted to flux densities with the solid angle subtended 
by the beam of the telescope via 



2k B T 



(5) 



Follow ing the algorithm described in Mac Low et al.| 
( 1991b| ), we convolve the resulting image with the beam 
width and add some noise according to the telescope 
parameters. The parameters for the Very Large Array 
( VL A^ are given in Table [l] and the parameters for the 
Atacama Large Millimeter Array (ALMA^are given in 
Table [2] We model the synthetic observation using an 
effective Gaussian beam rather than doing a full-fledged 
simulation of aperture synthesis. 

2.3. Generation of Dust Emission Maps 

For the wavelengths with A < 0.3 cm observable by 
ALMA, in addition to the free-free emission we must also 
take into account continuum emission by dust particles. 
We use RADMC-3E0 to generate dust emission maps 
as well as maps of combined free-free and dust emission 
from the simulation data. 

6 http : / /www. via. nrao.edu / astro / guides / vlas / current / node 1 1 . html 

7 http://www.eso.org/sci/facilities/alma/observing/speciflcations/ 

8 http:/ /www. mpia.de/homes/dullemon/radtrans/radmc- 
3d/index.html 
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TABLE 1 

VLA TELESCOPE PARAMETERS. 



VLA band 


wavelength (cm) 


beam width (arcsec) 


sensitivity (mJy) 


4 


400 


24 


160 


P 


90 


6.0 


4.0 


L 


20 


1.4 


0.061 


C 


6.0 


0.4 


0.058 


X 


3.6 


0.24 


0.049 


u 


2.0 


0.14 


1.0 


K 


1.3 


0.08 


0.11 


Q 


0.7 


0.05 


0.27 



Note. — The VLA can observe at eight bands with wavelengths from 400 cm to 0.7 cm. The beam width is given in terms of the FWHM 
and measured in arcsec, the sensitivity is measured in in mJy. The integration time is lOmin. 



TABLE 2 

ALMA TELESCOPE PARAMETERS. 



ALMA band 


wavelength (mm) 


beam width (arcsec) 


sensitivity (mJy) 


3 


3.1 


0.034 


0.019 


4 


2.1 


0.023 


0.022 


5 


1.6 


0.018 


0.411 


6 


1.25 


0.014 


0.044 


7 


0.95 


0.011 


0.079 


8 


0.7 


0.008 


0.272 


9 


0.45 


0.005 


0.411 



Note. — ALMA can observe at seven bands with wavelengths from 3.1 mm to 0.45 mm. The beam width is given in terms of the FWHM 
and measured in arcsec, the sensitivity is measured in mJy. The integration time is lOmin. 



RADMC-3D is an AMR-based radiative transfer pack- 
age for continuum and line radiative transfer. It h as an 



interface for PARAMESH ( |MacNeice et al.||2QQQ[ ), the 
AMR grid library of FLASH. We use RADMC-3D for 
two tasks. First, to compute the dust temperature self- 
consistently, using the standard Monte Carlo method of 
Bjorkman fc Wood| ( |2001| ), combi ned with Lu cy's method 
of treating optically thin regions ( Lucy|1999 ). Second, to 
compute the images of the free- free and dust continuum 
emission by using it as a volume-rendering ray-tracer 
tool. RADMC-3D is the succe ssor of the RADMC code 
( Dullemond fc Dominik| [2QQ4) which has been used in 



numerous papers. 

RADMC-3D has been tested against the earlier 2D 
version of RADMC for various ID and 2D test cases. 
A detailed discussion of RADMC-3D will be published 
separately, but since this is the first scientific use of 
the code, we show the results of a simple test case 
here. It involves a simple ID spherically symmetric en- 
velope around a star. The density of the envelope is 
Pdust(0 = Po(j/lAlJ) -2 , where po takes the values 10 -15 
gem -3 for test case 1 and 10 -14 gem -3 for test case 
2. The inner radius lies at 5 AU, the outer radius at 
100 AU. The star has solar parameters, but we treat the 
stellar spectrum as a blackbody of T = 5780 K. For the 
opacity we use silicate dust spheres of 0.1 /im, using op- 
tical constants of olivine from the Jena databasq^[ but 
we artificially set the scattering opacity to zero in order 
to be able to compare our results to the results from a 
simple ID variable eddington factor dust radiative trans- 
fer code called TRANSPHEREE3 With RADMC-3D we 
now compute the dust temperature using the Monte 
Carlo method. We do this in two ways. First, we use 
a spherical ID grid, similar to what we use for TRAN- 

9 http://www.astro.uni-jena.de/Laboratory/Database/jpdoc/ 

10 http:/ /www. mpia-hd.mpg.de/homes/dullemon/radtrans/ 



SPHERE. Secondly, we use a 3D Cartesian AMR-refined 
grid, where the refinement is done with the criterion that 
the cells with centers having radii r > 5 AU from the star 
should have a size Ax > 0.2r. This is a relatively coarse 
resolution, meant to test the effect of low resolution on 
the results. The outcome of this comparison is displayed 
in Figure [I] showing the excellent agreement between the 
temperature profiles. 

3. RESULTS 

Our simulations follow the gravitational collapse of the 
initial massive clump and lead to the formation of high- 
mass stars. The collapse leads to the formation of a 
massive rotationally flattened structure, which is gravi- 
tationally unstable and fragments. Only a single star is 
allowed to form in run A. It accretes 72M in 145 kyr. 
In run B, three high-mass stars with M > 10M Q form 
within 70 kyr, and become the dominant source of ion- 
izing radiation within the star cluster. The interaction 
of the ionizing radiation with the infalling accretion flow 
leads to multiple effects observable in both spatial and 
spectral diagnostics. 

3.1. Time Evolution 

The most striking property of the resulting H II regions 
is their extremely high variability in time and shape. In 
the online material of Paper I we presented movies of ra- 
dio continuum maps from different viewpoints. The ra- 
dio maps were generated for VLA parameters at a wave- 
length of A = 2 cm, using a beam with full width at half 
maximum of / /14 and a noise level of 10 _3 Jy. The as- 
sumed distance was 2.65 kpc. All the movies show the 
continuous build-up and destruction of UC H II regions. 
The timescale for changes of more than 5000 AU in size 
can be as short as 100 yr. 

This flickering is caused by the accretion flow in which 
the sources are embedded. When the protostar passes 
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Fig. 1. — Temperature profiles for two different core densities, 
po = 10 -15 gem -3 (left), po = 10 -14 gem -3 (right), as gener- 
ated by TRANSPHERE, RADMC-3D on a spherical ID grid, and 
RADMC-3D on a 3D AMR grid. All profiles agree with each other. 



through dense, gravitationally unstable filaments in the 
accretion flow, they absorb its ionizing radiation, so 
that the gas above the filament recombines and cools 
down. Since the gravitational instabilities cause the ac- 
cretion flow to be chaotic, the interplay between the radi- 
ation feedback and the infalling material results in highly 
stochastic ionization and recombination processes in the 
surrounding gas. 

This effect is demonstrated in Figure [2] It shows 
dramatic changes in the H II region around the most 
massive star in run B. Between t = 0.6592 Myr and 
t = 0.6595 Myr (within 300 yr), a region with a diameter 
of ~ 6000 AU suddenly recombines. Changes like this not 
only affect the physical size of the H II region, but they 
can also alter their morphology. From t = 0.6668 Myr 
to t = 0.6671 Myr (again within 300 yr), the morphol- 
ogy of the UC H II region surrounding the most mas- 
sive protostar changes from shell-like to core-halo be- 
cause of a large-scale recombination event that clears 
the rim of the shell. The shielding by the filaments also 
controls how ionizing radiation can escape perpendicu- 
lar to the disk. For example, this reverses the cometary 
H II region around the star between t = 0.6524 Myr and 
t = 0.6534 Myr (within 1000 yr). The three examples 
given in Figure [2] indicate that the morphology of ultra- 
compact H II regions around accreting massive protostars 
depends sensitively on accretion events close to the pro- 
tostar. 

The flickering observed in the simulations also resolves 



( |Wood fc Churchwell||1989| ). Since UC H II regions are 
not freely expanding bubbles of gas that monotonically 
increase in size, their diameter does not depend on their 
age. An extreme version of the discrepancy between pro- 
tostellar mass and size of the H II region occurs in run 
A, where the 70M Q protostar has almost no visible H II 
region. It is totally quenched by the strong accretion 
flow. 

While the source in run A never stops accreting, the 
most massive stars in run B finally stop growing when the 
gas reservoir around them is fully exploited. As the den- 
sity surrounding the most massive star then drops, the 
hot, ionized gas can finally expand monotonicially, even- 
tually overrunning the second most massive star as well. 
Once there is no more high-density gas in its neighbor- 
hood, the flickering around the most massive star stops. 
The ionized gas continues to expand, forming a compact 
H II region with larger size and fainter emission than the 
preceding ultracompact phase. We stress that it is not 
the ionizing radiation that stops the accretion flow by 
driving away the surrounding gas, rather it is the subsid- 
ing accretion flow that allows the ionized gas to expand. 
We call this process fragmentation-induced starvation, a 
more detailed discussion of which we presented in Paper 
I. 

3.2. Morphology 

The extended H II regions found in the simulations dis- 
play a large amount of substructure. Equation ([2| shows 
that the emission from free-free transitions scales with 
the square of the number density of free electrons, n e . 
This explains the emission peak close to the protostar, 
where very dense gas in the accretion flow gets partially 
ionized. However, not all emission peaks are associated 
with stars. 

Figure [3] shows some examples. The upper left panel 
(ai) shows an H II region with a shell-like structure. The 
shell clearly exhibits a peak on its rim that is several 
1000 AU away from any nearby star. The shell is created 
by dense shocks running through the H II region, which 
are replenished by material from the accretion flow. The 
emission of this dense gas is what creates the shell. An- 
other example is shown in the lower left panel (a2) of 
Figure [3] It shows a dense blob of gas that is externally 
irradiated by a massive star and creates a peak that ap- 
pears to indicate the position of a second star. Obviously, 
peaks in emission maps are not an ideal guide to the co- 
ordinates of stars. 

The aforementioned shocks contribute largely to the 
emission seen in the maps. The middle panels (bi) and 
(b2) in Figure [3] show an edge-on view of the rotation- 
ally flattened structure of the star cluster. The upper 
panel (bi) shows that the most massive star has created 
a cometary H II region. The lower panel (b2) displays 
the same region 200 yr later. The ionizing radiation has 
blown away gas from the accretion flow close to the pro- 
tostar. This shock runs away from the star and creates 
a filament of strong emission across the H II region. The 
right-hand plots (ci) and (02) in Figure [3] show the same 
region face-on. From this viewing angle, the shock shows 
up as shell-like structure. This demonstrates that the 
shell does not trace the edge where the ionizing radia- 
tion hits the accretion disk, but rather shocks generated 
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Fig. 2. — Changes in the H II region around the most massive protostar in run B. The images in the lower panels show the H II region at 
a slightly later time than the images in the upper panels. The left-hand panels (ai) and (a2) show the recombination of ionized gas with a 
diameter of ~ 6000 AU within 300 yr. The midway panels (bi) and (b2) demonstrate how the morphology of the H II region changes from 
shell-like to core-halo in 300 yr. The right-hand panels (ci) and (02) present the reversal of a cometary H II region in 1000 yr. The box size 
displayed is 0.122 pc. Black dots and circles indicate the position of all protostars with their accretion radii in the image. The mass of the 
central protostar is indicated in the images. All maps are simulated VLA observations at 2 cm with an assumed distance to the star cluster 
of2.65kpc. 



from inflowing gas. The shell-like structures around ac- 
creting protostars can be interpreted as indirect evidence 
for the accretion process. 

The origin of the shell morphology changes when accre- 
tion ceases. Figure [4] shows the late-stage evolution of the 
star cluster. The most massive star has stopped accret- 
ing, allowing its H II region to begin to expand quickly 
into the ambient gas. The left-hand plots and show the 
expanding H II region face-on in the upper panel (ai) 
and edge-on in the lower panel (a2). Here, the strong 
shell-like emission clearly comes from the dense gas in 
the rotationally flattened structure around the protostar 
rather than from a shock launched by the protostar. 

While the accretion onto the most massive star has 
stopped and cannot directly affect the structure of the 
growing H II region any more, it can still be influenced by 
other stars that interact with the gas. Two such events 
occur in run B and are shown in Figure [4j The middle 
panels (bi) and (b2) show a time sequence of a 7.8M 
star approaching the rim of the shell. Its ionizing radia- 
tion is strong enough to create sufficient thermal pressure 
to blow away the rim of the shell from its direct neigh- 
borhood. The right-hand panels (03) and (04) show the 
same star 8200 yr later, when it has already entered the 
compact H II region. Now its ionizing radiation can freely 
expand. The gravitational attraction of the star is strong 
enough to pull along a dense stream of gas. This stream 
allows the star to grow in mass although it has entered 
the H II region filled with underdense gas. This suggests 
that deformed shells could be indicative of stars inside 
large-scale H II regions. 



3.3. Dependence of Morphologies on Viewing Angles 

The simulation data presented here offer the unique 
opportunity to look at the same H II region from dif- 
ferent viewing angles. We have already seen that the 
observed morphology depends crucially on the position 
of the observer. We investigate this observation in detail 
for some H II regions appearing in run B. 

Figure [5] displays an H II region around the most mas- 
sive star in run B. The first panel (ai) shows the re- 
gion face-on, and the successive panels show the region 
rotated around an axis in the plane of the rotationally 
flattened structure by 18° increments until a rotation of 
90° is reached. Hence, the last panel (ag) shows the re- 
gion edge-on. This sequence of images demonstrates that 
a region with a shell-like morphology from one viewing 
angle can have a cometary morphology from a different 
viewing angle. 

The transition angle at which the shell-like morphol- 
ogy turns into a cometary morphology is about 72° in 
this particular case. The online material for this article 
contains a movie of the volume-rendered (logarithmic) 
absorption coefficient of this region along with the cor- 
responding synthetic VLA image as the line of sight is 
rotating, illustrating the complex structure of the ion- 
ized gas in the H II region, and how it determines the 
morphology, depending on the viewing angle. 

Because the ultracompact H II region shape is asym- 
metric, we get a different result if we perform the rota- 
tion around a different axis. Figure [6] starts with the last 
frame of Figure [5] and successively rotates the re- 
gion by 90° around the polar axis. This means that the 
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Fig. 3. — Emitting structures around the most massive protostar in run B. The left-hand panels (ai) and (3,2) give two examples of 
emission peaks that are not directly related to stars. This calls into question the usual understanding that coordinates of stars should be 
directly related to emission peaks. The middle panels (bi) and (b2) show an edge-on view of the star cluster with a cometary H II region 
around the most massive star. Within 200 yr, a shock is launched at the protostar, runs through the H II region and creates a filamentary 
structure. The right-hand panels (ci) and (02) show the situation face-on. Here, the shock looks like a shell expanding from the protostar. 
This indicates that the appearance of H II regions is closely related to accretion events close to the protostar. Black dots and circles indicate 
the position of all protostars with their accretion radii in the image. The mass of the central protostar is indicated in the images. All maps 
are simulated VLA observations at 2 cm with an assumed distance to the star cluster of 2.65 kpc. 



view is edge-on for all times. At an angle around 36°, 
the cometary region develops strong shell-like structures. 
We already know that these filaments have their origin in 
gas blown away by the protostar. This region is of shell- 
like type since it is bounded by a dense filament. This 
means that shell-like regions can occur both at face-on 
and at edge-on view. In principle, the same holds true 
for cometary regions. They can also be observed face-on 
when the ionizing radiation is shielded anisotropically. 

We have repeated the same exercise for different H II 
regions and find that the transition angle between shell- 
like and cometary morphologies shows no systematic be- 
havior. Figure [7] shows two more rotations from face- 
on to edge-on. We show only the first and last im- 
ages as well as the transition angle. The upper panels 
(ai) to (as) show the transformation of a shell-like mor- 
phology into a cometary morphology at t = 0.6864 Myr, 
the lower panels (bi) to (03) show a similar transition 
at t = 0.6925 Myr. The central star has a mass of 
22.532M© and 23.025M o , respectively. The transition 
angle is about 36° in the former and 54° in the latter 
case. They depend on the details of the structure. A 
universal angle above which transition occurs does not 
exist. 

3.4. Morphology Statistics 
3.4.1. Types 

Morphologies of UC H H reg ions were cla s sified by 
Wood fc Churchwell| fll989| ) and |Kurtz et al] fll994| ) as 
shell, co metary, core-ha l o, sph erical, irregular and unre- 
solved. De Pree et al. ([2005) abandoned the core-halo 



morphology and introduced a new bipolar category for 
elongated H II regions. They argued for abandoning the 
core-halo type because most H II regions are surrounded 
by faint emission, which produces a halo around any H II 
region. Though this may be true, we find it useful to keep 
this morphological type for regions with a pronounced 
central peak and a fainter envelope. The presence of a 
pronounced envelope clearly distinguishes this morphol- 
ogy from the spherical type, which we also find. 

In addition, we do not require shell-like regions to be 
void of central peaks. Although the larger, late-time 
shells do indeed not have central peaks, the UC H II 
regions associated with accreting protostars do have cen- 
tral peaks because they ionize their own accretion flow. 
In fact, observations with high sensitivity and resolution 
do find centrally p eaked shells that we re previously clas- 
sified as spherical ( |Carral et al.||2002| ). We predict that 
more regions of this type will be found as observations 
with better resolution and sensitivity also become avail- 
able for more distant massive star forming regions. 

Figure [8] demonstrates the importance of resolution 
and sensitivity in identifying the correct morphology. 
The figure shows synthetic maps of the same shell-like 
H II region for VLA parameters at 2 cm, placed at dif- 
ferent distances to the observer. The shell disappears 
between 6 and 10 kpc, where both the spatial resolution 
and the noise level make it impossible to distinguish be- 
tween parts of the H II region and pure noise. On the 
other hand, the number of centrally peaked shells seen in 
our simulations might be reduced by including additional 
feedback processes like line-driven stellar winds or mag- 
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Fig. 4. — Compact H II region created by the most massive protostar in run B. The left-hand panels show this region face-on in the 
upper panel (ai) and edge-on in the lower panel (a2). The emission in the shell has its origin in the rotationally flattened structure around 
the protostar and not in a shock running through the H II region. The middle and right-hand panels (bi) to (04) show a time sequence 
of a second star interacting with the dense gas that bounds the region. In the middle panels (bi) and (b2), this star blows away a dense 
filament by its own ionizing radiation. The right-hand panels (ba) and (b4) demonstrate what happens when it enters the low density 
region. Its gravitational field pulls a dense stream of gas behind it. Broken-up shells could thus be a helpful observational signature to 
locate stars inside H II regions. Black dots and circles indicate the position of all protostars with their accretion radii in the image. The 
mass of the central protostar is indicated in the images. All maps are simulated VLA observations at 2 cm with an assumed distance to 
the star cluster of 2.65 kpc. 
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Fig. 5. — H il region in run B around a star with 22.956M0 at t = 0.6907 Myr. The view in the first panel (ai) is face-on. The subsequent 
panels show rotations by 18° increments around an axis in the plane of the rotationally flattened structure. The last panel (a,e) shows 
the region edge-on. At an angle of about 72°, the morphology has turned from shell-like into cometary. All maps are simulated VLA 
observations at 2 cm with an assumed distance to the star cluster of 2.65 kpc. 
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Fig. 6. — H il region in run B around a star with 22.956M0 at t = 0.6907 Myr. The first panel (ai) is identical with the last panel (ae) 
of Figure [5] The region is successively rotated around the polar axis, so that the view is edge-on for all angles. The morphology changes 
from cometary to shell-like at about 36°. All maps are simulated VLA observations at 2 cm with an assumed distance to the star cluster 
of 2.65 kpc. 
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Fig. 7. — H il regions in run B that show different transition angles from shell-like to cometary. The upper panels (ai) to (as) show an 
H II region at t = 0.6864 Myr around a 22.532M star, the lower panels (bi) to (b 3 ) at t = 0.6925 Myr when the star has 23.O25M . The 
polar transition angles are 36° and 54°, respectively. All maps are simulated VLA observations at 2 cm with an assumed distance to the 
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netically driven outflows, which might be able to create a 
cavity around the protostar. If these processes would be 
strong enough to thin out the dense accretion flow suffi- 
ciently to remove the peaks remains an open q uestion. 
The new bipolar type is not well defined. 



TABLE 3 

Percentage Frequency Distribution of Morphologies 



De Pree 



et al. (2005) gave only one example for this new cat 
egor y, wher e the bipolar shape is not very distinctive. 
Churchwell (2002) required an hourglass shape for the 
bipolar morphology, which is not present in their exam- 
ple. Although we do find morphologies that look bipolar 
in the simulations, there are only very few of them, and 
their features are not very pronounced. The small num- 
ber of bipolar reg i ons is in agreement with observations 
( |Churchwell|[2002| |De Pree et al|[2005| . All of the bipo- 
lar regions could equally well tall into one of the other 
classes, which is why we do not take this category into 
account. 

The sensitive dependence on viewing angle and the 
high time variability of the H II region caused by the ac- 
cretion flow is sufficient to produce morphologies of any 
described type in a single simulation. Figure [9] shows 
maps from run A with only one ionizing source. The 
displayed shell-like and core-halo morpholgies are face- 
on views, whereas the comet ary H II region is viewed 
edge-on. As discussed above, the same morphologies can 
equally well be obtained at different viewing angles. It 
is apparent that the size of the H II region does not scale 
with the mass of the protostar. On the contrary, the ir- 
regular region corresponds to the largest protostar, but 
it is among the smallest H II regions. 

3.4.2. Comparison to Observations 

For comparison with th e ob servational survey s by 
Wood fc Churchwell| ( |1989| ) and |Kurtz et al.] ( [1994D , we 
perform a census ol morphologies in run A and run B. 
We select 25 snapshots from each simulation and view 
each one from 20 randomly chosen angles. This gives a 
total of 500 images per simulation. The different viewing 
angles take into account the fact that because of the ax- 
isymmetric geometry of the initial conditions, some mor- 
phologies preferentially occur at different orientations. 
For example, shell-like morphologies are found mostly 
face-on, and cometary morphologies edge-on. The set 
of different viewing angles, uniformly distributed on the 
unit sphere around the center of the computational do- 
main, avoids statistical biases by this effect. By the same 
token, the distribution of morphologies also changes with 
time. Since many stars in run B reach a mass of 10 M 
by the end of the simulation, the later times contain more 
spherical and unresolved H II regions than the beginning 
of the simulation. Since we do not know the geome- 
try and the evolutionary stage of the UC H II regions 
in the surveys, we assume randomly distributed orienta- 
tions and ages and thus average over different viewing 
angles and simulation snapshots to get a representative 
sample. We note that this somewhat overstates the con- 
tribution from later stages in our simulations, since with 
a standard IMF there are probably far more star form- 
ing events that terminate with the formation of low-mass 
OB stars than higher-mass ones. To achieve consistency 
with the observational analysis techniques we use contour 
plots to identify the morphological classes and follow the 
definitions given in |Wood fc Churchwell (1989). 

To guarantee that the viewing angles are distributed 



Type 


WC89 


K94 


Run A 


Run B 


Spherical/ Unresolved 


43 


55 


19 


60 ± 5 


Cometary 


20 


16 


7 


10 ± 5 


Core-halo 


16 


9 


15 


4 ± 2 


Shell-like 


4 


1 


3 


5 ± 1 


Irregular 


17 


19 


57 


21 ± 5 



Note. — Percentage relative frequencies of UC H n region mor- 
phologies in surveys and simulations. The table shows the mor- 
phology sta tistics of UC H II r egions in the surveys of |Wood fc| 
|Churchwell| |l989} (WC89) and |Kurtz et al.] fl994| (K94) as well 
as irom a random evolutionary sample rrom run A and run B of 
500 images for each simulation. For run B, the mean value and 
standard deviation of the relative frequencies from independent 
evaluations by the first four co-authors are given. The statistics 
for run A, in which only one massive star forms, disagree with the 
observations. 



evenly on the sphere, we use a rejection method to calcu- 
late the random angles ( |Press et al.|1986| ). We start with 
three random numbers x, z with a uniform distribu- 
tion on the interval [—1,1]. We then calculate the radius 
r = \/ x 1 + y 2 + z 2 and drop all points with r > 1. The 
remaining points are projected onto the sphere, X = x/r, 
Y = y/r and Z = z/r. The Cartesian coordinates 
(X, V, Z) of the point on the unit sphere are then con- 
verted into spherical coordinates resulting in the 
desired distribution of viewing angles. 

The results of our statistical analysis are presented in 
Table [U The numbers from run B are taken from evalua- 
tions performed independently by the first four authors, 
so that we can give mean values and standard deviations 
for the relative frequencies. The standard deviations re- 
mained small in absolute terms, never exceeding five per- 
centage points, although some are large in relative terms, 
particularly for the rarer types. 

Given the variation amongst the observational surveys, 
the results from the multiple sink simulation run B are 
quite consistent with the observational numbers. In par- 
ticular, we find that roughly half of the sample repre- 
sents spherical or unresolved UC H II regions, in agree- 
ment with the classical lifetime problem. Conversely, the 
results show that the relative numbers of spherical and 
unresolved H II regions in run A disagree by more than 
20 percentage points with the observational findings. 

This result clearly disagrees with theoretical models in 
which massive stars form alone. Since the protostar in 
run A grows very quickly, it cannot generate such a large 
number of strongly confined H II regions. Instead, lots 
of irregular H II regions form. The only way to get a 
large number of spherical and unresolved H II regions is 
the formation of a stellar cluster. This again shows that 
run B is a much more realistic model for massive star 
formation. 

3.5. Emission and Optical Depth 

It is interesting to study the appearance of H II regions 
at different wavelengths. Since the optical depth of the 
free-free radiation calculated from the simulation data 
is exactly known, we can investigate how morphological 
features of H II regions depend on the optical depth. As 
an example, we examine in detail the H II region from the 
upper panels (ai) to (as) of Figure [7] at t = 0.6864 Myr 
around a 22.532M star in run B. 
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Fig. 8. — Shell-like H il region of run B with a central peak placed at different distances to the observer. The synthetic 2 cm VLA maps 
show that the shell-like feature disappears between 6 and 10 kpc. To detect it at such large distances, observations with higher spatial 
resolution and sensitivity are required. Black dots and circles indicate the position of all protostars with their accretion radii in the image. 
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Fig. 9. — Different morphologies observed in run A, which has only a single ionizing source (central black dot). All morphologies found 
in surveys are present, depending on the simulation time and the viewing angle. The lower right panel demonstrates the problem with the 
bipolar category. The region is clearly elongated, but it also shows a shell-like structure. Since elongation alone seems to be insufficient to 
define a category on its own, we do not consider bipolar regions as a separate category. All maps are simulated VLA observations at 2 cm 
with an assumed distance to the star cluster of 2.65 kpc. 
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Figure [lO| shows synthetic VLA observations for wave- 
lengths from 20 cm to 0.7 cm assuming a distance of 
2.65 kpc (see Tab. [I]). The beam width decreases with 
wavelength, so that the H II region is increasingly bet- 
ter resolved. The corresponding images of optical depth 
are shown in Figure [TT] In agreement with the expected 
behavior evident from Equation ([2|, the optical depth 
decreases with decreasing wavelength. The broad peak 
at 20 cm appears not only because of the large beam, but 
also because the whole H II region is optically thick. At 
6 cm, only the right-hand side of the H II region is opti- 
cally thick, as is clearly visible in the VLA observation. 
At 3.6 cm, only a small region around the protostar and 
in an arc-like feature is optically thick, which leads to a 
shell-like H II region. The VLA maps at 3.6 and 2 cm 
look similar, but at the smaller wavelength the shell is 
not optically thick anymore. The one-to-one correspon- 
dence between emission features and optically thick re- 
gions breaks at 2 cm. For 1.3 and 0.7 cm, the optically 
thick region around the protostar disappears. The shell 
is still visible at these wavelengths. 

3.6. Millimeter and Submillimeter Maps 

At millimeter and submillimeter wavelengths one gen- 
erally expects to observe dust emission along with the 
free-free emission. This is not the case for the rather low 
mass core in the model presented here. Since the total 
initial gas mass in our model is only 1000 M , there is 
less than 10 M of dust in the whole domain during the 
simulation, under the assumption of a canonical gas-to- 
dust ratio of 100. This is not enough to produce any 
observable emission, because the free-free radiation is 
much stronger. Hence, the simulated ALMA maps of 
this model show only free-free emission, but not dust. 

We can estimate the maximum flux due to dust emis- 
sion as follows. Let Md us t be the total dust mass, d the 
distance to the source, n v the dust opacity and B V (T) 
the black-body spectrum for the temperature T. Then 
the total flux from all dust is approximately 



F v = ^M dust n v B v (T). 



(6) 



To arrive at an estimate for the upper limit of F V1 we 
set d = 2.65 kpc, M dust = 10M o , T = 1000 K and 
k v = 8 x lO^cnrV 1 at v = 187 GHz (A 1.6mm). 
This yields F v « 25.5 Jy or, assuming an UC H II re- 
gion with 0.03 pc diameter, 1.8mJy/beam for an effec- 
tive beam width of 0.018 arcsec. This is only a factor 
of 4 larger than the noise level for a lOmin integra- 
tion time of this ALMA band (see Table which is 
0.41 mJy. Though this calculation demonstrates that the 
small amount of dust in our simulation would already be 
very faint when free-free emission was absent, what really 
makes the dust emission invisible is its weakness com- 



pared to free-free emission (see also Section 3.7), which 
is independent of the telescope sensitivity. 

We caution, however, that the above calculation as- 
sumes a homogeneous distribution of dust across the H II 
region. This treatment neglects dust emission from an 
accretion disk around the massive star, which we cannot 
resolve in our simulation and may produce significant 
dust emission. In our numerical example, ALMA would 
be able to resolve these disks, if they exist, and detect 
their dust emission. A similar argument holds for disks 



around low-mass stars that are being exposed to ionizing 
radiation from a neighboring high- mass star, which may 
life long enough to be o bserved ( Hollenbach et al.| |l994 



Richling fc Yorke||1998[ ). This caveat also applies to the 
spectral ene rgy d istributions in the submillimeter regime 
(see Section |3?7| ) . 

We demonstrate the absence of dust emission in Fig- 
ure [12] with the H II region shown in Figure [TSJ The 
bright shell that was visible at VLA wavelengtns com- 
pletely disappears as the wavelength gets shorter. In- 
stead, bright spots directly at the location of the massive 
stars appear. This emission is caused by partially ion- 
ized, dense gas of the accretion flow around these massive 
stars. None of the images shows any appreciable dust 
emission. 

3.7. Spectral Energy Distributions 

A useful diagnostic of UC H II regions is their inte- 
grated SED. For a medium at constant temperature, the 
integral in Equation Q for the brightness temperature 
leads to 

T(7v)=T(l-e-n (7) 

In the frequency regime where the emission is dominated 
by free- free, Equations ([2| and (|5| show that the flux 
F v of a homogeneous source typically grows as v 2 in the 
optically thick regime (small frequencies) and falls off 
like u~ 0A in the optically thin regime (large frequencies). 
For typical emission measures of EM = 10 8 pc cm -6 , the 
transition frequency for which r = 1 is v t ~ 5.3 GHz. 
This turnover frequency depends only weakly on the 
emission measure, v t oc EM ' 48 . 

However, many observed SEDs of H 1 1 regions do not 
behave in thi s simple man ner but show anomalous scaling 
exponents ( |Lizano|[2008[ ). In particular, there are SEDs 
that grow linearly with v over a frequency interval that 
can be as large as the whole VLA range (see Table m). 
Such anomalous scaling exponents can be re produced 



H H region models with ioni zed grad ient s (|Panagia fc 
[FgTTi1[T975l |Omon||1975l JFranco et al.||2000| [Avalos et al. 
| 2006[ |Keto et al.||2008| ), hierarchical clumps ( |lgnace fc 
Churchwell 2004) and additional dust emission (Rudolph 
' ™' 1 — - " 1 2004]^ 



et al.|199UMPratap et al.||1992| |Beuther et al.||2(JU4| |Keto 



et al. 2008). 

The data from our numerical simulations is certainly 
more realistic than the simple ionized gradient or hier- 
archical clump models. Therefore it is reassuring that 
we can reproduce the typical shapes of observed SEDs. 
We show two examples in Figure 13 The left-hand SED 
shows the full coverage from VLAover ALMA to IRAS 
frequencies for the H II region of Figures 10 and T2| The 
SEDs show the typic al shape of those reported in |Wood| 



Churchwell| ( |1989[ ). As already noted above, the free- 
tree radiation is by far dominant up to ALMA frequen- 
cies. For IRAS frequencies, however, dust emission is the 
dominant process. The right-hand SED shows another 
example that more nicely illustrates the typical scaling 
behavior expected for the free-free emission. The dotted 
lines grow oc z^ 2 , the dashed lines fall off ex ^~ 01 . 

One example of an SED with an anomalous scaling 



oc v (solid line) is shown in Figure 14 The result can 



be a combined effect of density gradients in the ionized 
gas as well as shadowing by clumps, just as in the sim- 
pler analytical models. We do not need any dust to pro- 



12 



z(J cm 




(ai) 


6 cm 




(a 2 ) 


3.6 cm 




( a 3j 




40.00- 


n 






® 






♦ 










d 


30.00- 






® 

® 


® 

® 




® 

® 


i 

® 




® 

® 


® 

® 




ijy/bear 






















20.00- 






2 cm 




(a 4 ) 


1.3cm 




( a s) 


0.7cm 




(a 6 ) 


ssion in 












®^ 


® 










emi; 


10.00- 






• 

box size 0.122 pc 




® 

® 


® 

® 




® 

® 


® 

® 






o.oo- 







Fig. 10. — H il region in run B around a 22.532M0 star at t = 0.6864 Myr for different VLA bands. The beam width decreases with 
decreasing wavelength, so that the H II region is inceasingly better resolved. The assumed distance to the observer is 2.65 kpc. Black dots 
and circles indicate the position of all protostars with their accretion radii in the image. 
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Fig. 11. — The same H II region as in Figure [To] The images show the optical depth of free- free radition at the different VLA bands. 
With decreasing wavelength, the H II region becomes dominantly optically thin. The color scale is chosen with a sharp break at optical 
depth unity, so that the transition between the optically thick and thin regimes can be easily identified. Black dots and circles indicate the 
position of all protostars with their accretion radii in the image. 
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Fig. 



12. — The H II region from Figure [To] at ALMA wavelengths, 
in Figure [To] At high frequencies, only bright spots close to the massive stars are visible. 



Note that the color scale is logarithmic, different from the linear scale 
These images show no dust emission since it is 
too weak compared to free-free emission. This is because our whole simulation box contains less than 10 Mq of dust. Black dots and circles 
indicate the position of all protostars with their accretion radii in the image. 
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Fig. 13. — SEDs from H II regions in run B. The plots show full 
SEDs from VLA to IRAS frequencies. The dotted lines scale oc v 2 
and the dashed lines scale oc v~ QA . The agreement of the free- 
free emission with the expected power laws is different for the two 
regions. Dust emission is visible only in the IRAS range. Missing 
crosses representing dust emission are below the flux range shown 
in the image. The assumed distance for both SEDs is 2.65 kpc. 

duce anomalous exponents, and in fact in our particular 
model, the dust emission is so low that it is unable to 
produce such gradients. Hence, our simulations repro- 
duce both the general shape of H II regions SEDs as well 
as the observed anomalous scaling behavior entirely with 
the computed distribution of ionized gas. 

4. CONCLUSIONS 

We have presented synthetic continuum observations 
of the H II regions formed in our collapse simulations 
of massive star formation described in Paper I. We find 
that the H II regions are highly variable in size and shape 
as long as the massive star continues accreting. Dense 
filaments in the gravitationally unstable accretion flow 
irregularly shield ionizing radiation from the central star. 
This shielding effect leads to a large-scale (~ 5000 AU) 
flickering of the H II regions on timescales as short as 
~ 100 yr. As long as the flickering continues, there is no 
direct relation between the age of the star and the size of 
the H II region. This result appears able to resolve the 
UC H II lifetime problem. 

Furthermore, we have identified the structures inside 
the H II regions that produce the continuum emission. 
We find that emission peaks are not necesarily related 
to the positions of stars. Instead, strong emission is pro- 
duced by partially ionized gas from the accretion flow 
that enters the H II region. These filaments can be ion- 
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Fig. 14. — Anomalous SED from an H II region in run B. The 
dotted line shows a scaling oc v 1 , the dashed line a scaling oc i/ -01 
and the solid line a scaling oc v 1 . The free-free emission of this 
region grows anomalously with a spectral index around unity across 
the full VLA and most of the ALMA coverage. This scaling is 
produced solely by density inhomogeneities, not by dust emission. 
Missing crosses representing dust emission are below the flux range 
shown in the image. The assumed distance is 2.65 kpc. 

ized and blown away by thermal pressure. Since the con- 
tinuum emission is determined solely by the structure of 
the flow field around the massive star, the appearance of 
the H II region depends strongly on the angle from which 
it is viewed. For example, we find that the same H II re- 
gion can be classified as shell-like or cometary, depending 
on the position of the observer. 

We have evaluated the distribution of apparent H II re- 
gion morphologies. The multiple sink simulation run B 
reproduces the observed relative frequencies reported 
from surveys, including the high fraction of spherical or 
unresolved regions. The single sink simulation run A, 
however, fails to reproduce the morphology statistics. 
This is because the single star grows so quickly that it 
is statistically unable to produce a significant fraction 
of the smallest H II regions. The formation of a whole 
stellar cluster is necessary to get the relative frequencies 
of the smallest H II regions right. This analysis provides 
strong evidence against models of massive star formation 
where all or most high-mass stars form in isolation. 

In addition to the H II region morphologies, we can also 
reproduce the different SEDs characteristic of UC H II re- 
gions. We find that our initial gas mass of 1000 M is 
too small to produce observable dust emission at VLA 
or ALMA wavelengths. Nevertheless, ultracompact H II 
regions at different times and places in our models show 
SEDs with both regular transitions from optically thick 
to optically thin spectral slopes as well as anomalous 
scaling with a spectral slope around unity. These slopes 
are entirely produced by inhomogeneities in the density 
structure, with no contribution from dust. However, 
follow-up simulations with more massive initial clumps 
will be required to completely pin down the role of dust 
in observed SEDs. 
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